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Abstract. GYOTO, a general relativistic ray-tracing code, is presented. It aims at 
computing images of astronomical bodies in the vicinity of compact objects, as well 
as trajectories of massive bodies in relativistic environments. This code is capable of 
integrating the null and timelike geodesic equations not only in the Kerr metric, but 
also in any metric computed numerically within the 3+1 formalism of general relativity. 
Simulated images and spectra have been computed for a variety of astronomical targets, 
such as a moving star or a toroidal accretion structure. The underlying code is open 
source and freely available. It is user-friendly, quickly handled and very modular so 
that extensions are easy to integrate. Custom analytical metrics and astronomical 
targets can be implemented in C++ plug-in extensions independent from the main 
code. 
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1. Introduction 

The first developments of general relativistic ray-tracing date back to the 70s, with 
works regarding the appearance of a star orbiting around a Kerr black hole pQ, the 
derivation of an accretion disk's emitted spectrum in terms of a transfer function [2J, 
and the computation of the image of an accretion disk around a Schwarzschild black 
hole [3]. 

Since this period, many ray-tracing codes appeared in the literature, mainly 
dedicated to the computation of spectra and particularly of the 6.4 keV iron line 
[U El El El El EE] or to the computation of images of accretion structures [HI [T31 El [TO] . 
Some authors are primarily interested in simulating phenomena of quasi-periodic 
oscillations in the vicinity of black holes [131 E], the trajectory and appearance of stars 
orbiting around compact objects [T5J, [16] , the magnetohydrodynamical effects occuring 
in the surroundings of the black hole [T7] [T8]. the collapsar scenario of long gamma-ray 
bursts and the impact of neutrino pair annihilation [19, 20J, or the polarization of the 
propagating radiation [2T1 [22j [23] . In this context, one may wonder about the relevance 
of a new ray-tracing algorithm. There are two main reasons to do so. 

Firstly, it is important when using an algorithm to have access to its source code in 
order to handle it properly, and if needed to develop it according to one's needs. In this 
perspective, the code presented in this paper, GYOTO (General relativitY Orbit Tracer 
of Observatoire de Paris), will be available freely (see section [1]). Moreover, as GYOTO 
is written in C++, it is very easy to add new structures to the existing code, by using the 
object oriented aspect of the language. Those new structures, representing additional 
metrics or astronomical objects, can be compiled as plug-in extensions, independently 
from the Gyoto code itself. Although many ray-tracing codes are used in the literature, 
they are rarely made public (a recent exception is the Geokerr code [10]; note also 
the public release of GeodesicViewer for computing and analyzing individual geodesies 

mm)- 

Secondly, as compared with all other existing ray-tracing codes, GYOTO has a 
specific feature: it is capable of integrating geodesies in non-analytic, numerically 
computed metrics. This allows GYOTO to handle objects much more diversified 
than the standard Kerr black holes. In particular, the possibility to handle non-Kerr, 
numerical metrics allows to test the no-hair theorem of general relativity as investigated 
by [2S] in the context of analytic non-Kerr metrics. This kind of study with non-analytic 
metrics has not been done so far and will be at hand thanks to GYOTO. 

Basically, GYOTO consists in launching null geodesies from an observer's screen, 
that are integrated backward in time to reach an astrophysical object emitting radiation. 
Once the photon gets inside the emitting object, the equation of radiative transfer is 
integrated along the computed geodesic in order to determine the value of the emitted 
specific intensity that will reach the observer. 
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2. Ray-tracing in the Kerr metric 

2.1. Integration of geodesies 

The first goal of GYOTO is to integrate geodesies in the vicinity of rotating black holes. 
In this section, we will thus consider ray-tracing in the Kerr metric. 

2.1.1. Geodesic equations We use the standard Boyer-Lindquist (x a ) = (t,r,8,(p) 
coordinates to express the metric according to [27] : 

ds 2 = g a p dx a dx^ 

2Mr\ , 2 AMar sin 2 9 , , £ , 2 „ 1/l2 
dt 2 dt dp + — dr 2 + £ d# 2 



£ / £ A 

, 2 2 2 Ma 2 r sin#\ . 2 2 
+ r 2 + a 2 + sin 2 0d^ 2 , (1) 



where M is the black hole's mass, a its spin parameter (aM being the total angular 
momentum), £ = r 2 + a 2 cos 2 9 and A = r 2 — 2 M r + a 2 . 

For a particle of mass /i and 4-momentum p a , two constants of motion are provided 
by the spacetime symmetries: the energy "at infinity" E = —p t and the axial component 
of the angular momentum L = p 9 . The constant nature of E and L results respectively 
from the spacetime stationarity and the spacetime axisymmetry. In addition, the specific 
structure of Kerr spacetime gives birth to a third constant of motion, the Carter constant 
[28]: 

Q = p 2 9 + cos 2 9 [a 2 (fi 2 - E 2 ) + sin" 2 9 L 2 ] . (2) 

The problem is thus completely integrable by using the four constants (/i, E, L, Q). 

Following [TS], we use a Hamiltonian formulation for the equations of geodesies, 
which consists in using the variables (t,r,9,<p,pt,p r ,pe,Pip) to express the equations of 
motion according to 

r = —p r (3b) 
=^Po (3c) 

* = - ^1: (34) 

Pt =0 (3e) 
A\' 2 / 1 V 2 { R + AO 



^=-l2£j^-l2£j^ + l^A£-) ^ 
A\" 2 / 1 \ e 2 . f FL + AO 



pV = 0, (3h) 
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where the superscripts ' and 9 denote differentiation with respect to r and 9 respectively, 
a dot denoting differentiation with respect to proper time (for timelike geodesies) or to 
an affine parameter (for null geodesies). In addition, the following abbreviations have 
been used: 

9 = Q - cos 2 [a 2 (/i 2 — E 2 ) + sin" 2 6L 2 ] (4) 
R = [E(r 2 + a 2 ) - aL] 2 - A [//V + (L - aE) 2 + Q] . (5) 

2.1.2. Implementation Let us consider the integration of a single null geodesic by 
GYOTO. The initial conditions are the position of the observer and the direction of 
incidence of the photon. These quantities allow to determine the tangent vector to the 
photon's geodesic at the observer's position. 

The integration is then performed backward in time, from the observer towards 
a distant astrophysical objects. This backward integration is of particular use when 
the target object has a large angular size on the observer's sky: a large fraction of the 
computed geodesies will hit the object. If the integration was done from the object, 
most of the geodesies would not reach the observer, leading to a waste of computing 
time. However, owing to its modular nature, Gyoto could be fairly easily extended to 
support the forward ray-tracing paradigm. 

The equations of motion are solved by means of a Runge-Kutta algorithm of fourth 
order (RK4), with an adaptive step (see e.g. [29]). The integration goes on until one of 
the following stop conditions is fulfilled: 



• the photon reaches the emitting object (however, see section 2.2.1 for the particular 
treatment inside optically thin objects), 

• the photon escapes too far from the target object (what "far" means depends of 
the kind of object considered), 

• the photon approaches too closely to the event horizon. This limit is typically 
defined as when the photon's radial coordinate becomes only a few percent larger 
than the radial coordinate of the event horizon. However this limit must not be 
too large in order not to stop the integration of photons swirling close to the black 
hole before reaching the observer, thereby creating higher order images. 

In order to ensure a proper integration, the conservation of the constants of motion 
is continuously checked, and the solution of the integration is modified to impose 
this conservation. This choice of modifying the computed coordinates to ensure the 
conservation of constants is a consequence of our choice of Hamiltonian formulation for 
the geodesic equations. The usual form of Kerr geodesic equations (see e.g. [2TJ ) , which 
directly enforces the conservation of constants, involves square roots that lead to sign 
indeterminations. The Hamiltonian formulation allows to get rid of this problem [T5] . 
at the expense of slightly modifying the solution as explained below. 



The constancy of E and L is directly enforced by equations (3e) and (3h). The 



constancy of Q is imposed by modifying the value of 9 given by the RK4 algorithm, 
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according to ^ and (3c): 

= ±LJq_ cos 2 ^2(^2 _ £2) + sin -2 L 2] ? ( 6 ) 

where the ± sign is chosen according to the sign of 9 resulting from the RK4 algorithm. 
Finally, the constancy of the 4-momentum's norm is ensured by modifying the f 
coordinate furnished by the RK4 algorithm. For instance, when considering a null 
geodesic, the value of r is chosen to enforce the relation 

g tt i 2 + 2 g ttp i ip + g„ r 2 + g ee 9 2 + g w y? 2 = 0. (7) 

Thus 



± I 9tt i 2 + 2 gttp i ¥ + gee 2 + 5W ^ ^ 

9rr 

where the ± sign is chosen according to the sign of r resulting from the RK4 algorithm. 
However, in order not to depart too much from the RK4 output, the maximum 
modification of 9 and f is limited to 1% of their initial values. 

The numerical accuracy of GYOTO is investigated in the Appendix where a 
convergence test is given. 



2.2. Radiative transfer and spectra computation 

2.2.1. Radiative transfer Once the null geodesic, integrated backward in time from the 
observer's screen, reaches the emitting object, the equation of radiative transfer must 
be integrated along the part of the geodesic that lies inside the emitter. In order to 
perform this computation, two basic quantities have to be known at each integration 
step: the emission coefficient j v and the absorption coefficient a v in the frame of a given 
observer (typically the observer comoving with the emitting matter). These coefficients 
are related to the specific intensity's increment dl u as along the geodesic according to 

dl v = —a v l v ds and dl u =j l ,ds, (9) 

where ds is the element of length as measured by the observer. The impact of scattering 
is not taken into account by the code. However, it is possible, if needed, to include in 
the absorption and emission coefficients the effect of scattering. Moreover, the modular 
nature of GYOTO makes it possible to develop a specific scattering treatment that 
would be plugged into the existing code. 

As shown in [30], the relativistic equation of radiative transfer is 

— = S-AX, (10) 

where A is an affine parameter along the considered geodesic, and where the invariant 
specific intensity X, the invariant emission coefficient S and the invariant absorption 
coefficient A have been used: 
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v being the frequency of the radiation. These quantities are invariant in the sense that 
they do not depend on the reference frame in which they are evaluated. 

Let us now consider the reference frame comoving with the fluid emitting the 



radiation. The invariant equation ( 10 ) becomes in this frame 



^ — ^em [ji>em ^cm ^tm] > (12) 

where u em is the emitted frequency of the radiation. It is straightforward to obtain the 
following relation: 

dA v cm = ds em , (13) 

where ds em is the amount of proper length as measured by the emitter. This expression 
leads to: 

= j v - a v I v (14) 

J J "cm "em "cm \ / 

QS em 

where all quantities are computed in the emitter's frame: this is the standard form of 
the equation of radiative transfer in the emitter's frame. 



Equation (14) can be immediately integrated between some value sq where the 
specific intensity is vanishing (the "opposite" border of the emitting region) and some 
position s: 



I v (s) = J ™V{-J i a»(s"W') Us')ds'. (15) 
Provided that the absorption and emission coefficients are furnished at each position 



of spacetime, (15) allows one to compute the integrated specific intensity. This will be 



illustrated in section 12.31 

2.2.2. Spectra computation When considering an optically thick object, the 
computation of a spectrum is quite straightforward. Let us consider the spectrum 
emitted by a geometrically thin, optically thick accretion disk around a Schwarzschild 
black hole. Following [6], the specific intensity emitted at a given position of the disk is 
written 

Iu em oc <5(z/ cm - z^) e(r) (16) 

where 5 is the Dirac distribution and where e(r) follows a power law with parameter p: 

e(r) oc r- p . (17) 

By means of the invariant intensity X = J^/z/ 3 , the specific intensity observed by a 
distant observer can be related to the emitted specific intensity according to 

Iu ohs = 9 3 Iu cni , (18) 
where z/ Q b s is the observed frequency, and 

g = —. (19) 
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Figure 1. Emission line of a disk around a Schwarzschild black hole of mass M seen 
under an inclination of 45° . The internal radius of the disk is r = 6 M and the external 
radius r = 10 M. The power law parameter is p = —3. Only the contribution of the 
first order image is taken into account. This result can be compared with figure 3 
of 0. 



The observed flux F v is then related to the observed specific intensity according to 

d^ Dbs =/, obs cos0d^, (20) 

where Q is the solid angle under which the emitting element is seen, and 9 is the 
angle between the normal to the observer's screen and the direction of incidence. The 
GYOTO screen is a very simplistic model for a physical telescope equipped which a 
spectro-imaging camera. It this model, the screen is assumed to be point-like. A pixel 
on this screen corresponds to a direction on the sky, just like a pixel on a camera detector. 
Our screen object is defined, among others, by the field-of-view it covers on the sky, the 
number of samples (or pixels) to cover this field-of-view, and the spectral properties of 
the detector (number of spectral channels and their wavelength). The screen covers a 
solid angle Q on the sky and each pixel screen corresponds to a small solid angle around 
a given direction of incidence of the photons. Hence 

^oba = ^"obB.poel COs(0 pixe i) 6fy xe l, (21) 

pixels 

where I ! / oba)P ixei is the specific intensity reaching the given pixel, # P i xe i is the angle between 
the normal to the screen and the direction of incidence corresponding to this pixel and 
5Q pixe i is the element of solid angle covered by this pixel on the sky, defined as the total 
solid angle covered by the screen divided by the number of pixels: 

<5^pixcl = -r: — , (22) 

-< "pixels 

where / is the angle between the normal to the screen and the most external incident 
direction of photons on the screen. Figure [I] shows the resulting emitted spectrum for a 
Schwarzschild black hole seen under an inclination of 45°. 
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Figure 2. Image of a star at the ISCO of a Kerr black hole with spin parameter 
a = 0.9 M. Left: ray-traced image; the radial coordinate of the observer is r = 100 M. 
Right: some specific rays corresponding to the primary (red), secondary (green), 
tertiary (blue) and quaternary (magenta) images. The various black stars correspond 
to the source position at the emission of the photons and the black squares depict the 
black hole and the observer; the axes are labelled in units of M. 



When considering an optically thin object, the spectrum can also be computed 
quite easily. It is possible to compute maps of specific intensities for a given panel of N 
observed frequencies z/ obs j, 1 < i < N, by using (15), and to compute the observed flux 
at these frequencies by using (21). Illustrations of such computations will be given in 
section I 



2.3. Implemented astrophysical objects 

This section describes the various astrophysical objects that are currently implemented 
in GYOTO as potential targets for ray-tracing. 



2.3.1. Star GYOTO can compute the image of a moving star, orbiting around a Kerr 
black hole. The model for the star is very simple, only the timelike geodesic of the star's 
center is computed, and the star is defined as the points whose Euclidian distance to 
the center is less than a given radius R, i.e. the points whose Cartesian-like coordinates 
(x = r sin 9 cosip, y = r sin 9 simp, z = r cos 9) satisfy: 

{x - x c ) 2 + (y- y c ) 2 + {z- z c ) 2 < R 2 , (23) 

where (x c , y c , z c ) are the stellar center's coordinates. No internal physics nor tidal effects 
are taken into account. 

Figure [2] shows the resulting imag^f] for a star orbiting on the innermost stable 
circular orbit (ISCO) of a Kerr black hole of spin parameter a = 0.9 M. Four particular 
rays are depicted on the right side of the figure, in order to show the trajectory followed 

| Here, an image is defined as a map of specific intensity /„. 
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Figure 3. Image of a geometrically thin infinite accretion disk around a Schwarzschild 
black hole. 



by photons responsible for the primary, secondary, tertiary and quaternary images. The 
corresponding geodesies swirl around the black hole zero, one, two or three times before 
reaching the object. 

The capacity of GYOTO to compute the trajectories of stars around a Kerr black 
hole will allow to develop relativistic orbit-fitting scripts, in particular by using the 
Yorick package presented in section |4| 



2.3.2. Thin infinite disk A basic object implemented in GYOTO is the geometrically 
thin infinite accretion disk, with interior radius corresponding to the ISCO [31]. The 
expression of the emitted flux is given in [T2], and is proportionnal to the emitted specific 
intensity provided the emission is isotropic, yielding 



cx 



(p 2 - 3) P 5 



P 



V2 



In 



(3 - 2V2 



P 



+ V3 



P 



(24) 



where p = y/r/M, M being the central black hole's mass. The resulting image is 
depicted in figure [3] for the case of a Schwarzschild black hole and an inclination angle 
of 70°. It is in good agreement with the results of [T2] . 



2.3.3. Rossby wave instability in a thin disk In a geometrically thin disk, a Rossby 
wave instability (RWI) can be triggered if the density profile exhibits an extremum 
for some value of the radius trwi [32]. Density waves will be emitted away from the 
extremum region, with vortices appearing in the vicinity of r^wi- Such a disk has been 
implemented in GYOTO, the density profile being computed by means of the VAC 
code (see [33]). The observed flux emitted by such a disk is computed according to 
the prescription in |34J, as a function of the surface density and radius. The result is 
presented in figure El which can be compared with those obtained in |34j . 
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Figure 4. Image of a geometrically thin disk subject to a Rossby wave instability. 

2.3.4- I° n torus Ion tori are geometrically thick, optically thin accretion structures 
with constant specific angular momentum. They are the optically thin counterparts of 
the optically thick Polish doughnuts [35j EH] • 

For a torus made of a perfect fluid, it can be shown |35j that the energy-momentum 
conservation law along with the assumption of constant specific angular momentum leads 
to: 

^ = V M W (25) 
p + e 

where p is the pressure, e is the proper energy density of the fluid and the potential W 
is related to the covariant component u t of the fluid 4- velocity by W = — In \u t \. The 
isobaric surfaces thus coincide with the equipotential surfaces of W. 

The cross-sectional shape of the equipotential surfaces is given in figure 1 of [35] . 
It appears that one particular surface has a cusp at some r = r crit : it crosses itself in 
the equatorial plane. Equipotential surfaces contained inside this critical surface are not 
connected to the central object, thus matter cannot be accreted and swallowed by the 
black hole. It is thus assumed that the torus physical surface coincides with this critical 
surface. The central point, r = r central , coincides with the innermost equipotential 
surface and to the point of maximum pressure. 

It is convenient to introduce the dimensionless parameter 

w = . (26) 

Wccntral ~ W cr it 

The potential (and the pressure) increases continuously between r crit and r centra i. Thus, 
w varies from (cusp) to 1 (center) inside the torus. Outside, w can take any other 
value (positive or negative). 

The physics of the radiative transfer used for the ion torus is based on [37] and 
is presented in detail in [3B] ■ Figure [Bl shows the image of an ion torus around a Kerr 
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Figure 5. Image of an ion torus around a Kerr black hole of spin parameter a — 0.5 M 
with trivial radiative transfer (constant emission coefficient, no absorption). 



50 - 
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-50 50 
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Figure 6. Image of a synchrotron emitting ion torus around a Kerr black hole of mass 
M = 41O 6 M0, spin parameter a = 0.5 M and placed at the Galactic center, as seen 
by an observer on Earth. See also [38] , 

black hole with a spin parameter a = 0.5 M, with trivial radiative transfer included. 
This image can be compared with figure 5 of [39] . Figure [6] shows the image of such 
an ion torus at the Galactic center as seen by an observer on Earth, with radiative 
transfer including the emission of synchrotron radiation (see [3H] for the details of this 
computation). 

Figure [7] shows the synchrotron emitted spectrum of the previous ion torus. This is 
an example of the computation of the spectrum of an optically thin object by GYOTO. 
Astrophysically interesting perspectives of this kind of computation will be presented 
in EH]. 
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logio(v (Hz)) 

Figure 7. Synchrotron emitted spectrum by the ion torus of figure [6] for different 
spin values. The spin parameter is a = (solid black), a — 0.5 M (dashed blue) or 
a = 0.9 M (dotted red). See also [35]. 

3. Ray-tracing in a numerically computed metric 

3.1. Introduction 

Many spacetime configurations have been obtained via numerical relativity. Most of 
them have been computed within the 3+1 formalism of general relativity, which relies 
on the slicing of the four-dimensional spacetime by a family of (three-dimensional) 
hypersurfaces (see e.g. piOl HT1 W2\). The GYOTO code can perform ray-tracing in such 
a numerically given metric, allowing to consider astrophysical objects beyond the Kerr 
black hole, such as rotating neutron stars, binary neutron stars or black holes, stellar 
core collapse, etc. 

A spacetime (M,g a p) is given in a 3+1 form if the following data are available: 
(i) a family of spacelike hypersurfaces (Et)tgR forming a foliation of M. and (ii) a set 
(N, f5\ jij, Kij) of field^|]on each hypersurface S t such that iV is a strictly positive scalar 
field, called the lapse function, j3 l is a vector field on £ t , called the shift vector, 7^ is 
a positive definite metric on S 4 , called the 3-metric or induced metric and is a 
symmetric tensor field on E t , called the extrinsic curvature tensor (see e.g. [101 HD 112] ) • 
If (x l ) are regular coordinates on each E t , varying smoothly from one hypersurface to 
the next one, then (x a ) = (t, x l ) is a valid coordinate system for the whole spacetime 
(J\A,g a p). From the knowledge of (N, (3\ 7^), one can reconstruct the spacetime metric 
g a/ 3 according to 

gap dx a dx? = -N 2 dt 2 + lij (dx i + p i dt)(dx j + ftdt). (27) 

§ Following standard conventions, Latin indices span {1,2,3}, whereas Greek ones span {0, 1,2,3}. 
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3.2. GYOTO treatment of 3+1 metrics 

Assuming that the metric is provided in a 3+1 form, two techniques of geodesies 
integration have been implemented in GYOTO. 

Four- dimensional method: From the 3+1 data (A, (3 l , 7^), reconstruct the 4-metric g a p 



according to (27); then compute the Christoffel symbols 4 T a of g a $ with respect to 



the coordinates (x a ) and integrate the geodesic equation in the standard 4-dimensional 
form: 

d 2 X a . „ dA^dA^ 



dAs +r Virir=°- (28) 

where A is either the proper time (timelike geodesies) or some affine parameter (null 
geodesies), the geodesic being defined by x a = X a (\). 

Three-dimensional method: The geodesic is searched in the form 

x i = X^t), (29) 

i.e. the geodesic is parametrized by the coordinate time t and not by A. In [13], we have 
derived a system of differential equations for the functions X l (t), which involve only the 
3+1 quantities (N, /3\7y, K^j) and their spatial derivatives (denoted by <9j = d/dx l ) : 

dX i 

— = NV i - ft (30a) 
dV i 

= N[V l (V j dj In N - K jk V j V k ) + 2K i j V> - 3 T l jk V j V k ] 

(J. L 

-Y j djN - V j djft . (30b) 

Here the 3 r* jA , are the Christoffel symbols of the 3-metric 7^ with respect to the 
coordinates (x l ). The auxiliary variable V i (t) = N~ 1 (dX t /dt + is the linear 
momentum of the considered particle divided by its energy, both quantities being 
measured by the Eulerian observer (i.e. the observer whose worldline is normal to 
the Ef hyp ersurf aces) |43] . 

3.3. Discretization 



The differential equations (28) or (30) are integrated by means of a fourth-order Runge- 
Kutta algorithm. This requires the values of the 3+1 fields (N , j3 l , 7^ , Ky) at arbitrary 
points. If the numerical spacetime has been computed by means of a spectral method 
[H] . this does not cause any great difficulty since by essence such a method deals 
with fields and not with values at grid points (see below). If the numerical spacetime 
arises instead from a finite difference method [HI H2], the fields (Af, 7^- , Kij) are 
known only at the points of the grid used for the computation. One should then devise 
some interpolation procedure to get the values of the fields (A, (3 l ,7^, K^) and their 
derivatives at the points required by the Runge-Kutta algorithm. 
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Currently the geodesic integration is implemented in GYOTO only for numerical 
metrics computed by means of spectral methods. We are using spherical-type 
coordinates (x l ) = (r, 6, ip) on Et and divide the range of r in n computational domains 
(n > 2) : r G [0,r ], r G (r ,ri], r G (r n __ 2 , +oo). The spectral method consists in 
choosing finite numbers N r , Ng and N v of degrees of freedom in r, 9 and <p respectively 
(typically N r $^ ~ 20) and in expanding any scalar field / onto a polynomial basis: 

N v -l Ng-1 N, — 1 

f(t,r,e,p)= E E4 fe (t)%(e)e ifc (e)$ fe (^), (3i) 

fc=0 j=0 i=0 

where £ G [0, 1] (r < r ) or G [—1, 1] (r > r ) and is related to r by 

'r/ro for r < ro 

2r - r £ - r^_i 



re - r e -i 
_ 2r n _ 2 
r 



for r^_! < r < r#, i G [1, n — 2] 



(32) 



for r > r 



n-2 



and Ojfc and are (trigonometric) polynomial functions. Working with the 
Lorene library [35], we choose 

cos(mfcy?) for /c even 
sin(mfc(^) for k odd, 

[&/2] (integer part of k/2), 



*fc(<^) = 



%(0 



cos(j6») 


for mfc even 


sin(j6») 


for mk odd, 




for j even and r < ro 


T 2l+ i(0 


for j odd and r < ro 


r,(0 


for r > ro, 



(33) 
(34) 
(35) 



(36) 



Tj being the Chebyshev polynomial of order z, defined by the identity cos(ix) = Tj(cosx). 



Note that the (y9-basis defined by (33) is nothing but a Fourier basis, the field / being 



obviously periodic in <p. The choice (35) for the 6*-basis is motivated by regularity 
conditions associated with spherical coordinates (consider for example x = r sin 9 cos tp 
(m = 1) and z = rcosB (m = 0)). The choice (36) for the £-basis is motivated by 
the good properties of the Chebyshev polynomials with respect to the approximation of 



functions [33]- An alternative choice could have been Legendre polynomials. In (36), the 
well defined parity of the polynomial for r < r follows from the regularity conditions 
associated with spherical coordinates around r = 0. Finally note that thanks to the 
choice (32) for £, the last domain extends to r = +oo, which is reached for a finite value 
of £, namely £ = 1. 



GYOTO: a new general relativistic ray-tracing code 



15 



In view of (31), it is clear that within the spectral method, the scalar field / is 
entirely described by the N r x N$ x N v coefficients fijkif)- For instance, the radial 
derivative of / is evaluated as 

ar it ^V-lJVfl-l N r -1 

-£(t,r,0,<p) = ^EEE fiAt)Kj(0 Q AO)M<P), (37) 

fc=0 j=0 i=0 

where the derivative -R^(£) of the polynomial Rij(£) can be expanded as 

JV r -l 

= 5^ a iiP R P{3+^)- ( 38 ) 
p=0 

The index (j + 1) in the right-hand side reflects the change of parity induced by 
the derivative and known coefficients which depend on the actual family of 

polynomials. 

In the course of the Runge-Kunta procedure for the geodesic integration, when 



the value of one of the 3+1 fields is required at a given point, one uses formula (31) to 
evaluate it. If the numerical spacetime is stationary, the coefficients fijkif) do not depend 
upon t and are part of the known data. In a dynamical spacetime, one has to evaluate 
fijk{t) from the known data, which are the values fijk{tj) of the coefficients at a series 
tj of discrete times — the coordinate times at which the numerical spacetime has been 
computed. To this aim, an interpolation at order 3 is used, using the 4 neighbouring tj 
to compute the field value at t, by means of Neville's algorithm (see e.g. [29J). 



Once the coefficients fijk(t) are known, formula (31) provides a very accurate value 
of /, with an error decaying exponentially with the numbers of degrees of freedom N r , 
Nq and N v [44J. As a result, a number of degrees of freedom of order 20 is generally 
sufficient to reach the maximum relative accuracy allowed in double-precision computing 
(~ 1CT 14 ). Note that formula (31) involves a number of arithmetic operations which is 



of the order of N r X Ng x N 9 . Consequently the ray-tracing in numerical metrics is far 
more time consuming than that in analytical ones. 

3.4- An illustrative example 

To illustrate the computation of geodesies in a numerical metric, we consider a test-mass 
particle orbiting around a relativistic rotating fluid star. Contrary to the black hole 
case, the metric around the star is not known analytically (except when the star is not 
rotating, where it reduces to Schwarzschild metric). It has to be computed numerically. 
For this purpose, we have used the Lorene/nrotstar code [15], which is based on a 



spectral method in spherical coordinates (as described in section 3.3) and solves the 



Einstein equations according to the framework presented in (46]. The computation of 



the timelike geodesic is performed via both methods exposed in section 3.2 The results 



are compared in figure |8l revealing a very good agreement between the two methods. 
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Figure 8. Orbit of a massive particle around a rotating relativistic star, obtained as a 
timelike geodesic in a numerical metric resulting from the Lorene/nrotstar spectral 
code. The computation has been performed either by integrating the 4-dimensional 
geodesic equation p8| (solid red line) or by integrating the 3+1 geodesic system (30) 
(dotted blue line). The black circle is the central star's equatorial circumference. 
The axes are graduated in units of 10 km. The central star has a gravitational mass of 
2.1 Mq, is rotating at a frequency of 700 Hz, its equatorial coordinate radius is 15.2 km 
and the ratio of the polar to equatorial radii is 0.76. The orbiting massive particle is 
restricted to the equatorial plane, with its Z coordinate limited to \Z\ < 1 m. The 
radial coordinate of the star oscillates between a minimum value of around 38 km and 
a maximum of 210 km. 



4. GYOTO code in a nutshell 

GYOTO is documented in-depth on its homepage http://gyoto.obspm.fr. Here we 
simply summarize its most salient features. 

Several interfaces are provided to use the GYOTO code: 

• a command line utility (gyoto) allows ray-tracing a single frame. The scenery is 
specified through a XMljjj] file and the output frame is stored in the astronomy 
standard FITS^jjfile format; 

• an extension package (ygyoto) for the Yorickj^] interpreter allows using GYOTO 
interactively from a command prompt and writing complex scripts. A graphical 
user interface (gyotoy), based on this Yorick package, allows to interactively trace 
a single timelike geodesic in Kerr metric (see figure [9]); 

• all the core functionnalities of GYOTO are packaged in a C++ shared library 

|| http://www.w3.org/XML/ 

If http://fits.gsfc.nasa.gov/ 

+ http:/ /yorick. sourceforge.net/ 
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Figure 9. Screen capture of GYOTO's graphical user interface depicting an equatorial 
orbit in a quasi-extremal Kerr metric. It is possible to interactively change the 
parameters of the orbit thanks to the tabs on the right column. 

(libgyoto) and can therefore be easily accessed from the code. 

GYOTO has been designed with modularity in mind. A simple plug-in system 
allows users to easily add new metrics and new astrophysical objects without modifying 
or even recompiling GYOTO. The standard objects and metrics provided in GYOTO 
are also implemented in plug-ins, so examples are readily available. 

The GYOTO code itself is not (yet) parallel, but the ray-tracing concept is 
inherently parallelizable. Both the gyoto utility and the ygyoto Yorick extension 
package are designed to split image computation in regions and can be run from standard 
job submission systems such as torque, thereby performing an effective parallelization. 

GYOTO has been tested on recent versions of Linux and Mac OS X and should run 
fairly easily on any UNIX-like system. The code is available online at the following URL: 
http://gyoto.obspm.fr. A user guide is available there, that gives details about the 
code architecture, as well as about the software prerequisites needed to run GYOTO. 

5. Conclusion 

We have developed a new ray-tracing code, GYOTO. The code is public, and is 
designed in order to be easily handled by the end user, even if not specialist of ray- 
tracing. GYOTO is able to integrate null and timelike geodesies in the Kerr metric, 
as illustrated above for various astrophysical objects surrounding black holes, that are 
already implemented. Written in C++, GYOTO is particularly adapted to be developed 
according to the specific needs of the user, by virtue of its object oriented syntax. 

A specificity of GYOTO as compared with other ray-tracing codes is its ability to 
handle numerical metrics resulting from 3+1 numerical relativity computations. The 
development of this aspect of the code will certainly lead in the near future to interesting 
and new astrophysical results. 
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The near future of GYOTO will be devoted to the development of new astrophysical 
objects and to the diversification of the handled numerically computed metrics. The 
possibility to parallelize the code by using GPU will also be investigated. 
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Appendix: Convergence test 



Figure |A1| shows a convergence test of the GYOTO algorithm. We have considered 



a thin infinite accretion disk, as described in section 2.3.2 surrounding a black hole 



with different values of spin parameter. The observed flux was computed according to 



equation (21) for different values of the total number of pixels of the GYOTO screen. 



The flux reference value is defined as the flux obtained with a screen of iVpi x = 2000 x 2000 
pixels. The percentage of error, as compared to this reference value, is then determined 
for different values of N pix . 



Figure |A1| shows that GYOTO converges very quickly to very small errors. The 
error is already below 1 % for iV pix = 100 x 100. With N pix = 1000 x 1000, the value 
used for all the figures shown in this article, the error is less than 0.05%. 

The numerical accuracy of GYOTO is thus very satisfactory. Regarding the speed 
of GYOTO, the computing time necessary to obtain figure [3] is around 15 min on a 
laptop under Mac OS X on one core of a 2.4 GHz Intel Core 2 Duo processor. The 
image has iVpi x = 1000 x 1000, the observer being at a distance r = 100 M from the 
black hole. 
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parameter a = (black square), a = 0.5M (blue triangle) or a = 0.9Af (red diamond), 
as a function of the total number of pixels of the GYOTO screen. The percentage of 
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